cx = [-1, 0, 0, 1, 1, -1];
cy = [1, 1, 0, 0, -1, -1];

cx = cx;
cy = cy;

[x, y] = generateBoundaryNode(cx, cy, 30);

domain2D = createDomain(x, y);

dim = size(domain2D.x, 2);

b = assemblyBoundaryVector(domain2D);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

res = 30;
scale = 0.9;

[px1, py1] = meshgrid(0:res, 0:res);
[px2, py2] = meshgrid(0:res, 0:res);
[px3, py3] = meshgrid(0:res, 0:res);

px1 = reshape(px1, (res+1)*(res+1), 1) / res * scale - scale / 2;
px2 = reshape(px2, (res+1)*(res+1), 1) / res * scale - scale / 2;
px3 = reshape(px3, (res+1)*(res+1), 1) / res * scale - scale / 2;
py1 = reshape(py1, (res+1)*(res+1), 1) / res * scale - scale / 2;
py2 = reshape(py2, (res+1)*(res+1), 1) / res * scale - scale / 2;
py3 = reshape(py3, (res+1)*(res+1), 1) / res * scale - scale / 2;

shift = 0;

px1 = px1 - 0.5 + shift;
py1 = py1 + 0.5 + shift;
px2 = px2 - 0.5 + shift;
py2 = py2 - 0.5 + shift;
px3 = px3 + 0.5 + shift;
py3 = py3 - 0.5 + shift;

rx = [px1; px2; px3];
ry = [py1; py2; py3];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% [A, Z] = assemblyMatrix(domain2D);
% 
% theta = A \ b;
% 
% u = evaluateInteriorPoint(rx, ry, domain2D, theta);
% 
% F = evaluateInteriorDeformationGradient(rx, ry, domain2D, theta);
% 
plotStrain( domain2D, rx, ry )

% plotDeformation(domain2D, rx, ry, b(1:2:2*dim), b(2:2:2*dim), u(:,1), u(:,2));